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Abstract 

We analyse a p + ip-wave pairing BCS Hamiltonian, coupled to a sin- 
gle bosonic degree of freedom representing a molecular condensate, and 
investigate the nature of the BEC-BCS crossover for this system. For a 
suitable restriction on the coupling parameters, we show that the model is 
integrable and we derive the exact solution by the algebraic Bethe ansatz. 
In this manner we also obtain explicit formulae for correlation functions 
and compute these for several cases. We find that the crossover between 
the BEC state and the strong pairing p + ip phase is smooth for this model, 
with no intermediate quantum phase transition. 
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1 Introduction 



Progress in cold atom physics has yielded many studies into the nature of the 
BEC-BCS crossover [I]. Early theoretical accounts emphasized the need to 
study Hamiltonians which explicitly incorporate coupling between Cooper pairs 
of atoms and bosonic molecular modes [2]. Several works extended this approach 
to the case of p-wave paired systems [3] , a scenario that is experimentally acces- 
sible pi] . Currently there is substantial interest in p + ip-w&ve paired systems [5] , 
which has been primarily motivated by the seminal work of Read and Green [B] 
who illustrated the topological distinctions of the quantum phases occuring in this 
setting. Our objective here is to study a p + ip-w&ve pairing Hamitonian which 
is coupled to a bosonic molecular degree of freedom to investigate the BEC-BCS 
crossover in this context. Our approach is to employ exact Bethe ansatz methods 
for the analysis. 

There have been many exact analyses of the s-wave pairing reduced BCS 
Hamiltonian using the solution provided by Richardson [7J. These works were par- 
ticularly prevalent in the wake of experiments conducted on metallic nanograins 
[8]. A comprehensive understanding of the model's mathematical property of 
integrability has been developed [H] which has lead in particular to some in-depth 
investigations through the use of exact computation of correlation functions [10J. 
There have been efforts to extend these integrable methods to investigate models 
where there is coupling between Cooper pairs and bosonic molecular modes 
Generally, these examples fall into a class of generalised Dicke/Tavis-Cummings 
type integrable models |12j . They have the shortcoming that the pair-pair scatter- 
ing terms found in the Hamiltonians of [2] are not present, with only pair- molecule 
scattering terms appearing. 

More recently it has been established that an integrable model also exists for 
p + ip-w&ve pairing [T3HT6] . Integrability in this instance stems from a trigono- 
metric solution of the classical Yang-Baxter equation, in contrast to the rational 
solution associated with the integrable s-wave case. We will show below that an 
extension of this model through coupling to a bosonic degree of freedom, whilst 
maintaining pair-pair scattering interactions, is integrable for some restriction of 
the coupling parameter space. We will derive the exact solution of the Hamilto- 
nian's energy spectrum and certain correlation functions and use these results to 
study the BEC-BCS crossover. 

This paper is organized as follows. We begin Section 2 by introducing a 
general Hamiltonian describing ap + ip-w&ve pairing BCS model coupled to a 
bosonic molecular degree of freedom. Subsection 2.1 discusses the limiting case 
of the uncoupled system, in which the extreme limits of BEC and strong pairing 
BCS ground states are found. Subsection 2.2 establishes suitable constraints on 
the Hamiltonian's coupling parameters for which the system is integrable, while 
subsection 2.3 develops the exact solution via algebraic Bethe ansatz methods. 
The ground-state root structure of the Bethe ansatz equations is determined in 
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subsection 2.4, and based on these results it is shown in 2.5 that the ground-state 
wavefunction topology is trivial so no topological phase transition exists in the 
integrable case. Since the integrable case connects the extreme BEC and strong 
pairing BCS ground states, these belong to the same quantum phase. Section 
3 is devoted to the study of correlation function. Subsection 3.1 deals with 
one-point correlation functions and particular attention is given to the boson 
fraction expectation value. Subsection 3.2 deals with two-point functions and 
the boson-Cooper pair fluctuations are studied in some depth. Conclusions are 
summarised in Section 4. An Appendix on a mean-field treatment of the model 
is also included. 

2 Model Hamiltonian 

We consider a 2-dimensional p + ip-wave pairing BCS model coupled to a single 
bosonic degree of freedom where the Hamiltonian of the model is 

k 2 G 

H = 5b ] b + 7^ c i< c k - j Yl fa ~ ik v)( k 'x + ik' y ) c i c -k c * c -v 

k k^±k' 

- f J2 (fa ~ ik v) c l c -* h + h - c ) ■ C 1 ) 

k 

One sees that when 5 — K — 0, the Hamiltonian becomes the integrable p + ip 
pairing BCS model [13] with Ck and c k being destruction and creation operators 
of 2-dimensional polarised fermions, k and m the momentum and mass of the 
fermions and G a coupling constant which is positive for an attractive p + ip 
interaction. In the above Hamiltonian, the bosonic mode with destruction and 
creation operators b, W is associated to a zero-momentum molecular condensate. 
The interconversion between Cooper pairs and molecules is controlled by the 
coupling K. The sign of K is not important since it can be changed by the 
unitary transformation b — > —b. Included in the Hamiltonian is the detuning S 
which accounts for the energy splitting by a magnetic field due to the difference 
between the magnetic moment of the molecules and that of the Cooper pairs. 
Hereafter we set m = 1. This model is integrable if we set 5 = —F 2 G, K = FG 
with F being a free variable, which will be proved below. Before considering that, 
it is useful to first examine the ground-state phases of the uncoupled system. 

2.1 Limiting case of the uncoupled system 

Setting 5 = K = the Hamiltonian (JTJ, restricted to the Hilbert subspace where 
the bosonic degree of freedom is in the vacuum state, is the p + ip model. For 
the extended model ([I]), with 6 = K = on the full Hilbert space, the ground 



3 



state of the system is of the form 

10) = IVbcs) ® |iV 6 ) (2) 

where |?/>_bcs) is a ground state associated with the p + ip Hamiltonian and \N b ) 
is a bosonic number state. For the ground state we need to consider the optimal 
choice of the boson number N b which yields the lowest energy. Since the detuning 
is zero in this limit, the ground state will be one which provides the mimimum 
energy of \ipBcs) with respect to variations of the Cooper pair number. 

To elucidate the ground state structure in this limit we recall results from [T5] 
for the p + ip model, which has three ground-state phases called weak coupling, 
weak pairing, and strong pairing. Letting Nc denote the number of Cooper pairs, 
we set %c = Nc/ X as the filling fraction, and g = GC Throughout, 2C denotes 
the total number of momentum levels such that £ is the number of momentum 
pairs. The three phases are characterized by the constraints shown in Table 1 . In 
the weak coupling phase the ground-state energy is positive, on the Moore-Read 
line it is zero, and in all other cases it is negative. Ground states in the weak 
pairing and strong pairing phases, with filling fractions dual whenever 

x^ + x c = 1 — g^ 1 , with the two ground states having the same energy. The 
Read-Green state is self-dual. The Read-Green condition xq — (1 — g l )/2 gives 
the state with the lowest possible energy, for all g > 1, with respect to variations 
in xc- For g < 1, corresponding to the weak coupling phase, the lowest possible 
energy is given by the vacuum since all ground states with xc > have positive 
energy in this phase. The only phase for which the ground-state wavefunction is 
topologically non-trivial is the weak pairing phase [T5] . 



Phase 


Filling fraction xc 


weak coupling 


x c > 1 - g~ L 


Moore- Read line 


xc = 1 — g^ 1 


weak pairing 


(l-g- l )/2<x c <l-g- 1 


Read-Green line 


x c = (l-g- 1 )/2 


strong pairing 


x c <{l-g- 1 )/2 



Table 1.- Ground-state phases of the p + ip model. 



In view of the above we can determine the ground-state structure of ([I]) when 
5 — K — 0. We let x = x b + xc denote the filling fraction of the system, where 
Xb = Nb/C. If g < 1 all p + ip states with xc > have positive energy, so the 
ground state is obtained by choosing x b = x and xc = 0, giving a pure BEC state 
for (T5]) with zero energy. For g > 1 the p + ip ground states have negative energy. 
If x > (l — g~ 1 )/2 we choose xc = (1 —g~ 1 )/2 so the p+ip state is the Read-Green 
state, which has the minimum energy with respect to variations of xc- This then 
leaves x b — x — (1 — g^ 1 ) /2 so ([2]) is mixed. Finally if x < (1 — g~ v )/2 the p + ip 
state is in the strong pairing phase. The energy is miminised by choosing x b = 0. 
This leads to the classification shown in Table 2. 
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Phase 


9 


Filling fraction x 


\^BCS) 


N b 


BEC 


<?<! 


all 


vacuum 


N 


Mixed 


<7>1 


x > (l-g- 1 )^ 


Read-Green 


0<N b <N 


BCS 


9>l 


x<{l-g- l )/2 


strong pairing 






Table 2.- Ground-state phases of the Hamiltonian ([I]) for 6 — K — 0. 



To investigate the crossover between the BEC state and the BCS state we 
may start with g > 1 and K = 5 = in the Hamiltonian ([1]) so the ground 
state consists of the strong pairing p + ip state and the bosonic vacuum provided 
x < (l—g~ 1 )/2. By turning on K and 5 we obtain an interacting system of Cooper 
pairs and bosons. Next we vary g such that g < 1, and then turn off K and S. 
The ground state will now consist of the p + ip vacuum and a bosonic number 
state. The question we ask is whether the system experiences a phase transition 
as we pass from the strong pairing BCS state to the BEC state in this manner. 
Importantly, the coupling parameters can be varied such that the Hamiltonian 
remains integrable as we move between the BEC state and the strong pairing 
BCS state. 



2.2 Integrability conditions for the coupled system 

It is convenient to first perform a transformation on the Hamiltonian (OQ). We 
enumerate the complex momenta k = k x + ik y , with k y in the upper half-plane, 
by integers j = 1, ...,£. Implementing the canonical transformation 

kx %ky . . 

Sj CkC—k, Zj | "■ | ? 

we may rewrite the Hamiltonian (JT|) as 

H = 5N b + (1 + G)H - GQ ] Q - KQ% - KtfQ, (3) 

where we have defined 

N b = tfb, N J = s] Sj , j = 1,...,£, 

c c 

H = j2 z * N i> 3 t = E*i4 (4) 



provided we restrict to the subspace of the Hilbert space which excludes blocked 
states (see jS] for a discussion of the blocking effect). This restriction is sufficient 
to study the ground-state properties when the total fermion number is even. We 
use N to stand for the pair number operator which is the sum of the boson and 



pairing number operators; namely, N = N b + N c with N c = Ylj=i-^j- We note 
that iV commutes with the Hamiltonian (DD). This allows us to block diagonalise 
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the Hamiltonian into sectors labelled by the eigenvalues of N, which are non- 
negative integers. Hereafter we will adopt the practice to interchangably use the 
symbol iV to denote the pair number operator and its eigenvalues. 

Now we show that for a suitable restriction on the coupling parameters of (OQ) 
the model is integrable. The integrable manifold is defined by the relations 

5 = -F 2 G, K = FG (5) 

with F being a free variable. Under this constraint, the Hamiltonian ([1]) becomes 

H = -F 2 GN b + (1 + G)H Q - GQ ] Q - FGQ ] b - FGb^Q. (6) 

We will prove the integrability of the above Hamiltonian by using the Quantum 
Inverse Scattering Method [17]. Our approach is a generalisation of the method 
detailed in [T5] . 

Let V be the 2-dimensional f/ g (s/(2))-module and R(X) E End(V <g> V) the 
six-vertex solution of the Yang-Baxter equation 

acting on the three-fold space V <E> V <8> V. The i?-matrix, which depends on the 
spectral parameter A and the crossing parameter q, explicitly reads 

/ Xq 2 - \~ l q- 2 | \ 

A -A" 1 | q 2 -q- 2 



R(X) 



q 2 -q- 2 | A- A" 1 

\ j Xq 2 - \- l q- 2 J 



We construct the Yang-Baxter algebra by using the i?-matrix and the L-operator 
L(A) through the Yang-Baxter relation (YBR) 

R12 (A/ y) Lij (A) L 2 j (//) = L 2j { t i)L lj {X)R 12 {X/f Jl ). (7) 

Here L QJ (A) G End(V <g> V) is a 2 x 2 matrix of operators. In the framework 
of quantum integrable systems, the subscript a labels the auxiliary space, while 
entries of the matrix are operators acting on the j'th quantum space. 

A well-known L-operator is realised by the i?-matrix itself which, using local 
creation and destruction operators s, is expressed as 



L a i (A) 



\ q m-i) _ X -i q -m-V ( g 2 _ q - 



2 



Si 



(q 2 - q- 2 )s\ Ag-^- 1 ) - A _1 g (2JVi_1) 



where Ni is the local number operator with the definition Ni = s|sj. The op- 
erators s\, Si and are generators of the quantum algebra U q (sl(2)). In the 
2-dimensional representation they satisfy the relation 



Sj, s 



]]= 5^(1-2^). 
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A realisation of the L-operator using the g-boson algebra was given by Kundu [18J: 



\ q 2N b _ iA -l g -2(JV 6 +l) _ g fcr/4 ( ? 4 _ ? -4) l/2 & 
_ e -/4( g 4 _ g -4)l/2 6 f Ag -27V 6 + a -l g 2(W 6 +l) 



The subscript 6 in the L-operator L ab (X) stands for the bosonic quantum space. 
The local g-boson operators b q , b q and N b = b q b q have the following commutation 
relation 

? 2(2JV(,+1) + g -2(2iV b +l) 

[ 6 9> 6 J] = 2 i _ 2 • 

« ^ _|_ g 2 

It can be seen that when q — > 1, b q and fej become the usual bosonic destruction 
and creation operators b and b*. With the help of the mapping 

L ab {X) -e-/ 4 (g 4 - g-^VMiag (g 1 / 2 , g^ 2 ) ■ L ab (X) ■ diag (g 1 / 2 , g" 1 / 2 ) 

and the variable shift A — > — e~ i7r / 4 (g 4 — g _4 ) 1 ^ 2 A the L-operator, which still 
satisfies (ED, becomes 



Lo-b(A) 

where the elements are 



(Lcrbjll (L CT b)i2 
(L(jfe)21 (L ab ) 2 2 



to 



(A*)n = Ag^ +1 ) - (g 4 - g-*)X~ l g^ N ^\ 
(L ff6 )i2 = (g 4 - g^)b qi 
(L ab ) 21 = (g 4 - g- 4 )6j, 

(L CTfc ) 22 = Ag-^ +1 ) + (g 4 - g- 4 )A" 1 g^ +1 ). 

Now we define the monodromy matrix 

T CT (A) = g a L ab {Xz^)L aC {Xz c l )---L a2 {Xz^)L al {Xz^)) (8) 

with the diagonal matrix g a = diag(e~"*, e 4a )( cr ). Using the YBR ([7]), the following 
equation holds for the monodromy matrix 

R ap (X/fi)T a (X)T p (ii) = TMTrWRe^X/n). 

This relation ensures the commutation relation 

[t(A),t(/i)] = 0, VA,/i 

where t(X) is the transfer matrix defined by t(X) = tr a [T a (X)}. 
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Expanding the transfer matrix t(X) in orders of the spectral parameter A 

c+i 



i=-C-l 

we find that the coefficients commute with each other 

[t (l) ,t (j) ] =0 

for all In this manner we may construct an integrable system by using the 
coefficients ft\ The leading terms of the expansion are 

= ^flz-^j (e- ia q 2N - c+1 + h.c), 
t& = 0, 

= - (V f[ Z A E z ) (e-^q 2N - c+1 q^ N ^ + h.c.) 
\ i=i / j=i 



1 ) (e-* a q 2N - c+1 q-^ + V - h.c) 

i=i J 

+ (? 2 -?- 2 ) 2 U- i rR 1 



c 



i=l 

C / k-1 



x^ ¥ Je-V M+1 II q- m - 2) Sk*l + h.c. 

j<k \ l=j+l 



3=1 V 1=3+1 

Introducing the notation 

q = e ip ) /3 = rip, a-0(2N-C + l)=rit, 
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we define a Hamiltonian H by using the coefficient t i - c ^: 



i=l 

£ / , , , ^ x ., c 



sin2 M + W~2p \ 1 * 

sin 2 (2r/p) 3 \ 2 / 2sin 2 (2r/p) ^ 

+ 2^ 2 cot(2r/p) sin(r/(t + ApN b + 2p)) 



Ai-1 k-1 

„%>(4jVj-2) „ 

k a 3 



+ Y, z i z *[ e ~ int II e-^ Nl - 2 h k s) + e^ J] e^ 4 ^- 2 )4 S 

c / c 
+ 2cos(2 W ) z ^z^ e -^ +2 ^ +p ) J] e-^'-^ft^t 

j=l V Z=j+1 



_|_ e i»)(i+2pAr i ,+p) "Q e ir,(4Ar ; -2) & f s ^ j _ 



(9) 



Let G = 2p/t and F = 2zb- Taking the limit rj — > 0, we obtain the following 
Hamiltonian 



r?— >0 



2 sin (2r?p) ^ 7 ^-f J 16p 2 p 



£ £ £ 

= z ) N i ~ F2 ° Nb ~ G Y1 z i Zk + h - c ) ~ FG J2 z i + h - c ) ■ 

j=l j<k j=l 

(10) 

Utilizing (j3J) we find that ffTOj) is equivalent to fl6]). Therefore we have established 
that the constraint (J5} defines an integrable manifold in the coupling parameter 
space of (JTJ. 



2.3 Algebraic Bethe ansatz solution 

The eigenvalues of the Hamiltonian ffTUl) can be obtained by using the algebraic 
Bethe ansatz. Again, we follow the procedure of [15] and only present the main 
results. Rewriting the monodromy matrix T CT (A) (|H1) by using global quantum 
operators A(X), B(X),C(X) and D(X) defined by 



T(X) 



(MX) B(X) 
\C(X) D(X) 



the transfer matrix t(X) becomes 

t(X) =tr CT [T a (\)]=A(\) + D{\). 
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The Bethe states of the system are defined by 



N 



i=i 

where |0) is the vacuum state with the definition 

6|0> = Si \0) = 

for all i = 1, . . . , C. 

By using the standard algebraic Bethe ansatz method [TTJ , the eigenvalues of 
the Hamiltonian (Q are given by 



E = 24 cot(2,p) S in(„( t + 2p)) - ( ^iff^ ) £ 

2 sin 2 (2?7p) ^ sin(2r/p) ^ ^ ' 

Here, the parameters \ij (j = 1, 2, . . . , N) satisfy the following Bethe ansatz 
equations 



e 



N -1 _2 -1 2 

Taking the limit r) — >■ 0, we obtain the eigenvalues of ( fTOl ): 

JV 

E = {l + G)Y,ti 

3=1 



subject to the Bethe ansatz equations 

i A i, 2 — r 2 ~~ 



G- 1 + 2A^-£-l 4z 2 1 _A 2 



for j = 1, ... , N. For convenience, throughout the remainder of the paper we will 
simplify notation by making the substitutions 

A 2 H> A r //j i ^ //,-. 
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such that the Bethe ansatz equations take the form 

G- + 2iV- £ -l + 4j + f ^ = f-g_. (11) 

In the 77 — > limit the Bethe states and their dual states are defined by 

N 

l0({A})> = n C ( A .)l°>> (12) 

3=1 

N 

(0(M)i=(oin^)' 
3=1 

where C and B are the global creation and destruction operators given by 

C(A) -^ + ±^. ,13) 

j = l J 

B{») = 2 Jl + ^J^i_ (14) 

2.4 Ground-state root structure 

It is necessary to understand the character of the roots of ffTTj) which correspond 
to the ground state of the model. By adding an appropriate constant term to 
06]), all matrix elements are real and negative on each sector with fixed N. From 
numerical studies of ffTTT) we find solutions for which all the roots Xj are real and 
negative. From (IT51) . and an appropriate rescaling C(A) — > — C(A), we see that 
these roots give rise to an eigenvector with positive components. This eigenvector 
necessarily corresponds to the ground state as a result of the Perron-Frobenius 
theorem. This theorem also tells us that there is a unique solution set with the 
property that all Xj are real and negative. 

While we are unable to prove existence of a solution set with the property of all 
roots being real and negative in a general finite system, we can establish existence 
of such a set in the thermodyanamic limit. To analyze the thermodynamic limit 
of the model, £ — > 00, N — > 00 such that the filling fraction x = N/C remains 
finite, we follow the approach of [rSllTB] used to treat the strong pairing phase of 
the p + ip model. Making use of the following notations 

nr . 2z b G- 1 + 2N-C-1 
9 = GC, / = 7 =, g= 

and assuming that the ground-state roots fij become dense on an interval [a, b] 
of the negative real axis, the BAEs (ITT]) become the integral equation 

& *)_i4_p/V^=0 (15) 
e — ix n Li J a n — Li 
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where r(p) is the density of the roots, and p{e) is the density of the e = z 2 located 
on the positive real axis such that 



de p(e) = 1. 



The filling fraction x and intensive energy cq = lim E /jC are given by 



X 



dpr(p), e = / dppr(p). 



(16) 



Using standard techniques of complex analysis, the solution r(/i) of ([15]) is 



r(/i) 

R(p) 
S 



TV l 



de- 



S T 



LJO 



e — p)R(e) pj p? 



y/(ji- o)(//- 6), 



1 



2Vab 
P 



f 2 (a + b) 



Aab 



2Vab 



with the constraint 



q + 



2ab 



— Vab / de 



R(e) 



Evaluating ( TIBj) gives 
1 / 2 



a6 



e 



2 4Va6 



gp(g) 
fl(e)' 



2e: — a — b 
2R(e) 



(17) 



(18) 
(19) 



The value for e is obtained by solving equations (j!7|18p for a, b, and substituting 
these values into ffl9|) . Equations f|T7|) -f lT9|) are in agreement with mean-field 
results given in the Appendix. 



2.5 Topology of the ground- state wavefunction 

In the case of the p + ip model the ground-state phases as depicted in Table 
1 are independent of both the distribution of the momentum variables and the 
cut-off u, which is a consequence of the topological nature of the phases. For the 
analysis of ([1]) we consider that the momenta are fixed, so the parameter space of 
([1]) is three-dimensional with G, K, 5 as the variable coupling constants. For the 
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two-dimensional surface within the parameter space for which the Hamiltonian 
([!]) admits an exact Bethe ansatz solution, the ground-state roots are real and 
negative. We can use this property to show that the ground-state wavefunction 
is topologically trivial in the exactly solvable case. To this end we will adopt the 
winding number approach used in [13l[T5] for the p + ip model. 

The topological structure of a complex function <p(k) = <p x (\s) + i(p v (k) can 
be characterized by a winding number w. Consideration of the stereographic 
projection of the k-space domain of y?(k), and the stereographic projection of 
the image of </?(k), induces a map between Riemann spheres <p : S 2 —¥ S 2 . We 
adopt the convention for the stereographic projections that the point at infinity 
for both k-space and the image of y?(k) is associated with the north pole of the 
spheres. The winding number associated with tp is 

_ 1 f „ „ d kx p x d ky p y - d ky ip x d kx (p y 

W — — I UK X dKy y ; 7) ; ^7-9 • 

The key point to recognise is that a non-zero value of w can only occur if the 
north pole is in the image of (p, which is equivalent to the statement that y?(k) 
is divergent for some k. These concepts generalise to multivariate functions 
<p(ki, ...k A/ ). 

Now we turn to the ground state as given by f|T2|) and consider the expansion 

N 

\<p) = J2\^®\ Nb = N -ti 

j=0 

where each \ipj) is a state of j Cooper pairs expressible as 

hfc>= E ^( k i'-' k i) c L ct - kl - c k/-kjo}- 

ki,...,kj 

The possible pole structure of ^'(ki, k 3 -) can be deduced from the co-efficient 
terms of each C(fjLj), viz. 7(k) given by 7(k) = ik x — ik y )/ (/i — k 2 ). Since the fij 
are all real and negative for the ground-state these terms do not diverge for any 
k. This situation should be contrasted with the p + ip model where changes in 
the topology of the ground state occur exactly when some of the roots Hj vanish, 
in which case 7(k) diverges at k = p2J[T5]. Hence the functions ^(ki, •••,k 7 ) 
are topologically trivial. This leads to the conclusion that there is no topological 
phase transition in the exactly solvable case. As the exactly solvable case allows 
us to crossover from the strong pairing BCS state to the BEC state, these two 
states belong to the same topological phase of the ground-state phase diagram. 

3 Correlation functions 

Our conclusion that the crossover between the BCS and BEC ground states is 
smooth should manifest in the correlation functions of the Hamiltonian. Within 
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the framework of the algebraic Bethe ansatz, these may be computed exactly 
in terms of determinants of matrices whose entries are functions of the roots of 
(JTTT) . Following the calculations of [15] (see Appendix A. 2), based on results of 
Slavnov [19j, if the parameters fii satisfy the Bethe ansatz equations (1TTT) . the 
scalar products of states for arbitrary parameters Xj are 

S(M,{A}) = (<&({/,}) |$({A})} 

= uN 77 * -detgm,{\}) (20) 

with 

Q,m, {A}) = Rim - A,) ^ + E r,,.,^.,^ + £ 



Of particular interest is the case when ^ — Aj Vz, whereby 

5({A},{A}) = detG({A}) (22) 

with 



4^ 2 £ 7 2 - 2A 

G({A})„ = ^^^-^-^—^ 

Equipped with this result, we now proceed to calculate several forms of correlation 
functions. 

In general an m-point correlation function is defined by 

nw, 4, • • • , c {ad = m^wi ■ ■ ■ o({a})>, 

where e^. stand for any local pairing operators , s|. , or bosonic operators 
6, and iV?,, and the lower indices ij indicate the positions of the operators. With 
the help of the definition of the global creation and destruction operators (I13|ll4p . 
we may solve the inverse problem for the local operators through 

2 2 

a) = lim -C(v), sj = lim L B(y), 

&t = lim— C(v), b=lim—B(v). (23) 
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Below, instead of representing the local number operators Nb and N m in terms 
of global operators, we will use the following commutation relations to compute 
their correlation functions 

[N b ,C{\)] = j±b\ (24) 

Note that throughout we assume that the parameters {/ij} satisfy the Bethe 
ansatz equations (TTTj) . For the off-diagonal one-point functions below, the cardi- 
nality of the set {/ij} is one greater than that of {Xj}, while in all other instances 
they have equal cardinality. 

3.1 One-point correlation functions 

• The off-diagonal one-point function for the fermion pair creation operator 

^(M,4,{A}) = (0(M)I4I0({A})). 

Substituting the representation of ( )23|) into the above definition, we have 

lim v -^(<t>m)\c{\ 1 )...c{\ N ^)C{v)\Q). 



It is seen that the one-point function is a limit of the scalar product ( 1201) . Sub- 
stituting the scalar product into the above formula, we obtain 



pm, si, {A}) - n ^ - zl) det(r - (M ' {A})) 



with the elements of the N x N matrix given by 

CC(M, {A})) .. = (Q({p}, , j = l,...,iV-l, 



W(M,{A})) tj 



We also note that 

F(M,4,{A}) = F({A},s m ,M). 
• The off-diagonal one-point function for the boson creation operator is 

Fm,b\{\}) = (<j>m)\tf\<i>(W)) 
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Using a similar method as for the case of sj^, we obtain the one-point function 
as follows: 

F({fi}, b\ { A}) = N ^m^m — 

IL=i A * Uk<i ( x k - Xi) UkKii^i - Vk) 



with 



(%m, {A}))<, = {G{{ij}, {A})) y , j = l,...,iV-l, 

(r fe (M,{A}))^ 



2 • 



/' 

Similar to the case above we have 

Fm,b\{\}) = F({\},b,{ri). 

• To calculate the diagonal one-point function for the Cooper pair number 
operator N m , we consider only the functions 

F({A},iV m ,{A}) = (0({A})|iV m |0({A})). 

Using the commutation relation (|25p and the fact that iV m |0) = 0, we obtain 
F({A},iV m ,{A}) 

N 

= E rz^WWWAi) ■ ■ ■ C(Ai-i)4AA i+ i) . -.C(A^)|0> 

j=l ra 



(0({A})i0({A})}-(0({A})|0({A})} 

9 

v — zi 



+ E lim 2 \ — ^mmc(\i)...c(\^)c(v)c(\ J+1 )...c(\ N )\o} 



= detG({A}) - det (G({A» - Q m ({A})^ 
with the elements of the rank-one matrix Q m ({A}) being 

Q m ({A}) N 



z 2 



ij (Xi — Z 2 n ) 2 

In the above derivation, we have used the following property of determinants: 
If A is an arbitrary n x n matrix and B is a rank-one n x n matrix, then the 
determinant of A + B is given by 



det (.4 + B) = det A + detA(i) 



i=l 

where the elements of the matrix A.® are defined as 



A% = A aP for /? ^ i, 



16 



• Here we calculate the diagonal one-point function for the boson number 
operator N b 

F({A},iV f) ,{A}) = (0({A})|iV b |0({A})). 

Similar to the above case, using the commutation relation (|24"1) we obtain the 
one-point function for N b as 

N ?7 N 

F({\},N b ,{\}) = ^^(0({A})|6tJ]c(A fc )|O) 

j=l 3 k^=j 



A, 

= (0({A})|0({A}))-(0({A})|0({A})) 

N N 
j=l 3 k^j 

= detG({A}) - det (<5({A» - Q b ({\})) (26) 
with the elements of the rank-one matrix Q b {{\}) reading 

(4({A») .. = ^f. 

Through this last example we can compute the boson expectation defined by 

{ b) (0({A})|0({A})> • 
Substituting (1221) and (1261) into the above definition, we obtain 



det G({A})-Q 6 ({A}) 

W = 1 7^ \ '- (2T) 

det G({A}) 



and in turn the boson fraction expectation value (N b )/N which has previously 
been used to characterise BEC-BCS crossover properties [20] . 

Since the ground-state roots of the Bethe ansatz equations (jTTj) is the unique 
solution set which is real and negative, this makes for an efficient study of the 
ground-state features in finite systems. This is because the numerical solution of 
ffTTT) for negative real roots is very reliable, due to the uniqueness of a solution with 
this property. As an example, we take 2C = 900 momenta which arise in pairs k 
and — k. The distribution of the momenta is chosen as |k| = v2n, n = 1, ...,450, 
giving the cut-off as u = 30. Taking the total particle number as M = 300, 
corresponding to the filling fraction x = 1/3, we numerically solve for the ground- 
state roots of the BAEs f fTTj) to calculate fl27j) . In this sector the Hilbert space 
has dimension 1.96 x 10 123 . The results shown in FigJT] suggest smooth variation 
of the boson fraction expectation value for f,g>0, consistent with the absence 
of a phase transition. 
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£=450, N=150 



(Nb> 
N 




4 



Figure 1: The boson fraction expectation value (N^/N as a function of the 
coupling parameters /, g, as given by (127)) . The results shown are for a system of 
C = 450 momentum pair states and N = 150 pairs, giving the filling fraction as 
x = 1/3. The boson fraction expectation value shows smooth variation between 
the BCS ((N b )/N = 0) and BEC ((N b )/N = 1) extremes. 

3.2 Two-point correlation functions 

• We first determine the off-diagonal two-point correlation function for s\J) 
^({A},4&,{A}) = (0({A})|4&I0({A})). 
The commutation relation between operators b and C(A), viz. 

allows us to commute the bosonic operator with all C(Xj). Considering that 
b\0) = we obtain 

N 17 N 

j = l 3 fc-gLj 

= (0({A})|0({A})}-(0({A})|0({A})} 

_ 2 N 9 N 

+ lim ^Dr^l^IIWi ) 



Aj 



i=l fe^j 



detG({A}) - det (^({A}) - M m ({A})) (2* 
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with the elements of the rank-one matrix M m as 



V / %J Aj{Ai-z m ) 

The canonical two-point function with the definition 

^ (0({A})|0({A})) ' 

is expressible, using f !2S|) . as 



det G({A})-M n ({A}) 

(&4> = 1 * 7^ ^ (29) 

det \ G{{\}) 



Fig. [2] illustrates the behaviour of the ground-state two-point function (bs^) as a 
function of the coupling parameters /, g. In all instances there is a rapid decrease 
in the fluctuations as / — > 0, but they appear smooth nonetheless for f,g>0. 
Fig. [3] shows the scaling behaviour of these two-point functions. These results 
indicate that, for a fixed value of g, we may write 

(bsi)=MxJWn(f/C) 

where (p n (x, f) is finite and 9 n (f / C) has the property that 9 n (0) = 0. For the 
thermodynamic limit N, C — > oo with x = N/C we conclude that (bs^) — > 
0, which is one of the main assumptions underlying the mean-field treatment 
discussed in the Appendix. 

For the remainder of this subsection we calculate three more cases of two-point 
correlation functions. Although we will not numerically evaluate these examples, 
the formulae are included for completeness. 

• Here we calculate the two-point correlation function for s^Sn 
F({\},sls n ,{\}) = (0({A})|4 Sn |0({A})) 

Operating s n on the state |0({A})), we have 

N N 

Su\<P({X}n)) = s n Y[C(\p)\0) = Sn J] (Cn(X p ) + ^4) |0>, 

0=1 p=l 

where 

~ 2z b b ] •s—y zis l p _ z n 
UlA^j - — — + 2^T 3' A n~\ 75- 

A /3 A /3 _ Z l A ~ Z n 
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Figure 2: Ground-state two-point correlation function (bsl) (1291) . as a function 
of the coupling parameters / and g. These two pinjt functions represent the 
boson-Cooper pair quantum fluctuations in finite systems where (b) = (si) = 0, 
Here C = 100 and the distribution of the momenta is |k| = \/2n. The insets 
correspond to the cases (a) iV = 50, n = 1; (b) iV = 50, n = 100; (c) N = 150; 
n = 1; (d) N = 150, n = 100. For intermediate values 1 < n < 100 we have 
found that (&s£) has the same generic profile. It is apparent that increasing iV is 
associated with increasing fluctuations. 
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Figure 3: Relationships between the boson-Cooper pair quantum fluctuations 
(bsl) and the rescaled coupling parameter f/y/C with g = 100, for different 
pairs of (C,N). The momentum distribution has been rescaled as |k| = v/2n/ '£. 
Inset (a) shows the fluctuations for the case n = £ as functions of /. Here 
there are nine distinguishable curves corresponding to the choices of (£, N) where 
£ = 100, 200, 330 and N = 50, 100, 150. The remaining insets for (b) n = £, (c) 
n — £/2, (d) n = 1 illustrate that, with the inclusion of scaling factors, the pair 
quantum fluctuations can be expressed by functions of the variable f /y£. These 
cases also display data for the nine choices of (£,N), which is seen to fall on a 
single curve. The scaled fluctuations go to zero as / /VX — > 0. This indicates that 
the fluctuations vanish in the thermodynamic limit N, £ — > oo, with x = N/£ 
finite. 
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Bearing in mind that (s^) 2 = 0, s n |0) = 0, we find that the non-zero terms in 
the above relation combine to give 

TV N 

*n|0({A})> = $>^nCn(A,)|0> 

P=l ajtp 
N N 

= e^ii( c ( a «)-^)io) 

13=1 a^P 

N N N N N 

p=l a^p 13=1 a</3 l¥=<*,P 

With the aid of (15U|) . the two-point function reduces to 

N N 

F({X},sls n ,{\}) = ^^(0({A})|4II C ( A «)I ) 

0=1 a+f3 

N N N 

P=l a</3 1^a,P 

Now the two-point function for s^Sn has been simplified to sums of one-point 
functions of s^ m and two-point functions of s^s^. For the first term, we have 

N N N 

E A n(<P(W)\sl II C ( A «)I°> = E W(A/3 - z 2 m )] detf£({A}), 

0=1 /3=1 



where 



(^({A}))^. = (g({A})).., j?J3, 
f%({\}) 



ip (\ l -z 2 m ) 2 ' 



For the second term, we have 

N N N 

2 EE^(0{A})i44 n c ^\°) 

/3=1 a</3 ~t+afi 

, _ 2 w _ 2 \ N N N 

= lim lim Zm)[V Zn) J2J2 2A n A n^W)\C(u)C(v) J] C(A 7 )|0) 

P=l a<p 77^, P 

AT / TV 

= E [^"( A /3 - 4)] E^ det ^(W)) 

p=l V a</3 
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with 

«({A}))« = (G({A») , 
(S({A})) fa = tt-^. (^({A}))^ 



(A,-^) 2 ' V mnVL J "* (Ai-^) 2 ' 
_ a/3 = 2^ TO (A a - ^)(A /3 - zp 

We can therefore express this two-point function as 
F({A},4 Sn ,{A}) 

= E W( A * - £)] f det ^({A» - E / ^ det (^({ A }))) • 

/3=1 V a</3 J 

Writing the columns of the matrices and T^ n in vector notation, we have 

det(f^) = det^i,...,^-!,^,^!,...,^), 
det(f^) = det(G 1 ,...,G a _ 1 ^ n ,G a+1 ...G^ 1 ,W m ,Gp +1 ,...,G N y 

where Gj, j = 1, . . . , N denotes the N- dimensional vector with entries (Gj)i = G^ 
and W x (x = m, n) denotes the N- dimensional vector with entries 

(W x ) 



Note that we have refrained from detailing the explicit dependence on {A} in the 
above expression, and hope it is still clear to the reader. We focus our attention 
on simplifying the expression 



detT£ - K^detT* 

a</3 



mn 



for each permissible 8. Using properties of determinants, it is possible to establish 
that 

detT^ - £ KldetT^ = detX^ 

a</3 



where 



X 13 ) — G K jP — ?' < 8 

7 V \ A i - Z n) 



X 



Z 7 < 



fe) = Gy, J>0. 
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Therefore, the two-point function simplifies to 

N 

F({\},sls n ,{\}) = 

(3=1 

We remark that the above procedure for reducing the double sum of determinants 
to a single sum leads to a more compact expression compared to the analogous 
result in [T5] . 

• Two-point function of N m N n 

Now we consider the two-point function of N m N n 

N 

F({\},N m N n ,{\}) = (0({A})|iV m iV n J]C(A l )|O). 

i=i 

Commuting the number operators by using the commutation relation ( 1251) . we 
derive the correlation function as follows: 

N 

F({\},N m N n ,{\}) = (0({A})|AWV n n C ^)l°) 

i=l 

N N 

= E t^tWI^H^II^^I ) 

13=1 Afi m a^P 
N N N 

= ErTrErb^W II ^10) 

P =l A ? Z m a+fi A « Z n liz(x3 
N N 

f3=l oc±{3 
N N 

/3=1 a</3 

AT V N 

= E [ det ^ - det ^] + E E( J - - 4 a Jdet(f^), 

/3=2 (9=1 a</3 



.(A/j 



detXL({A}) 



where 

ja/3 _ z m z n{^a ~ Z m )(Aff ~ gn) 

mn " (^-^)(A a -A^) • 

At this stage it is worth pointing out the obvious fact that in the second term in 
the last line of calculation above, the summation never sees the (3 = 1 term, so 
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we can proceed by writing 

N 

F({\},N m N n ,{\}) = 

13=2 



N 

det(f£) + Y,(J% ~ ^n)det(f^) 

a<f3 



N 

-E det (^)- 

0=2 



For each (3, using similar techniques as before, the expression in square brackets 
above can be simplified to a single determinant, namely 

N 

det(7& + £(Jl - J&)det(7*£) = det(F^), 
where the matrix elements are 

( ^mn ) . . = + — Jmn)j\ 7<f\2' ^ ^ ^' 

\ = Z m 

k™)* ^ (\-zir 

Yin).. = G ij: j>/3. 
Once again using familiar properties of the determinant, we may also simplify 

N 

£det(^)=det(4n), 

0=2 

where the matrix A m has elements given by 

(A m )ii = Gn, 

(A m )ij = G i:j - G ij+1 , 1 < j < N, 
(A m ) m = j^j- 

In the above we have supressed the explicit dependency on {A} in each of the ex- 
pressions, as it should be clear. Therefore the two-point function can be expressed 
as a sum of N determinants 

N 

F({\}, N m N n , {A}) = det(Fl({A})) - det(A m ({A})). 

13=2 

• Two-point function of N b N m 
To compute the two-point function 

N 

F({\},N b N m ,{\}) = (0({A})|iV b iV m J]C(A,)|0), 

i=i 
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we need to use commutation relations both (1!M]) and (I25I) . The result is given by 

N N 9 AT 

F({A},iV b iV m ,{A}) = ^-^^-^(0({A})|4& t II C ^)l°) 

0=1 A P Zm a ^0 Aa j^a,/3 
N 2 N N 

=.^E^|SE{WW)R*) n <wi<» 

m ,8=1 p m a 7^«,/3 



iV 

' mn) 



} j } j ■Jmn ( ^ e ^(^mi 
0=1 a^0 
N 

EE( J -- J -)det(5^) 

0=1 a<0 



where 



ja/3 _ 2z b A )3 (A Q , - z m ) 

(Dm n )ij — Gij (j 7^ (^mn)ia ~ ~T2~' {^rnn)i0 



Using similar techniques as before, we can reduce the above to a sum of iV 
determinants. Doing this leads to 

N 

F({\}, N b N m , {A}) = det(ZL({A})) - det(A m ({A})) 

0=2 

where the elements of Z^ n are given by 

Z mn ^j — Gij + 2(Jj^ — Jmn)~^2i 3 ^ 

Z 



mn J -a ~ f X r2 X2 ' 
Zmn) .. = Gij, j > (3. 



4 Conclusion 

We have introduced a model that couples a p -Hp- wave pairing BCS Hamiltonian 
to a bosonic degree of freedom, and studied its properties regarding the BEC- 
BCS crossover. For a restriction on the coupling parameters, the model was 
shown to be integrable and the exact solution was derived by the algebraic Bethe 
ansatz. We found that the ground-state roots of the Bethe ansatz equations have 
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the property that they are the unique solution set such that all roots are real 
and negative. Using this result we reasoned that the ground-state wavefunction 
is topologically trivial, so the BEC-BCS crossover is smooth. This conclusion 
was supported by a study of the boson fraction expectation value, which was 
computed exactly in the Bethe ansatz framework. We also formulated expressions 
for a range of two-point correlation functions and used one particular example 
to study the boson-Cooper pair fluctuations. The range of correlation function 
expressions we have obtained provide ample scope for further studies along the 
lines of pig. 
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5 Appendix - Mean-field theory 

In the mean-field approach products of operators A and B are approximated by 



where the notation (■) stands for the expectation value. This approximation 
assumes that quantum fluctuations may be neglected. Formally, we define the 
fluctuations to be 



where A = 2G(Q) + 2K(b) is referred to as the gap. Since (152"]) does not commute 
with iV the Lagrange multiplier u, the chemical potential, has been introduced 
in order to tune the expectation value (N). 

*We omit the term GHq which becomes negligible in the thermodynamic limit which is 
discussed in Subsection 12.41 



AB « A(B) + (A)B - (A)(B) 



(31) 




A(GQ f + Ktf) 



(32) 
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We take the following form for the mean-field variational ground state 

= | a) <8> ® 1^2) ® •••• ® 1^/:), 

where |a) is the coherent state such that b\a) = a\a) with a = \a\e^ 6 C, and 
1?/^) are local states related to the pairing operators 

with |0j) the vacuum state of the momentum pair space labelled by j. Minimising 
the ground-state energy expectation value and imposing self-consistency of other 
operator expectation values leads to the following set of equations determining 
the values for a, the gap A, and the chemical potential v: 

KA . . 

a = 2(G5-Gv + K*y (33) 

r £ 2 

5 - V Z 



G5-Gis + K 2 



E / (34) 

, _ ,,\2 1 ,21 AI2 



i= i A /(zJ-z,) 2 + ^|A| 



2N-C- — — -— + 



2(G5 — Gv + K 2 ) 2 G5-Gu + K 2 
c 

3 1 1 3 1 



E / 1 ( 35 ) 

< * I/O \ O OlAlO 



J=1 w^-i/)a + 2?|A|> 



The ground state energy assumes the form 



1 * „ / 2z 2 + \A\ 2 -2u , 



2 



and we also find 



\Uj\ 




z)-v 





— V 


) 2 + z) 


\A\ 2 




z 2 


— V 






— V 


)> + z 2 


|A| 2 



< = n \ l-- F =J== \=(N j ). 



Let v 2 = ab and 2v — |A| 2 = a + b. It can be verified that taking the 
continuum limit for flMl) - fl36|) with the substitution (jSj) reproduces equations (fTTj) - 
ffT9|) . Moreover when (J5J holds, (|35|) - (|55)) lead to the result 



(N b ) \a\ 2 f 2 \A\ 2 + [" d£ p(e)(e + Vab) 



N N lQxu 2x 2x J R(e) 
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This quantity is a smooth function for /, g > 0, which can be shown by using the 
fact that a, b < 0. It displays discontinuous behaviour (in the first derivative) 
only in the non-interacting limit / — > when g = 1 or g^ 1 = 1 — 2x, which is due 
to level crossing. These correspond to the transition points of Table 2, and are 
visible in Fig. 1. The smoothness of the boson fraction expectation value for the 
interacting system in the thermodynamic limit is consistent with the absence of 
a phase transition between the BCS and BEC extremes. 
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